


00 
(N 



> 



Simulation of Cavity Flow 
by the Lattice Boltzmann Method 



Shuling Hou^'^, Qisu Zou^, Shiyi Chen^, Gary D. Doolen^, Allen C. Cogley^ 



C . 1 

j^ , Center for Nonlinear Studies and Theoretical Division, Los Alamos National Laboratory, 



Los Alamos, NM 87545 
2 



Department of Mathematics, Kansas State University, Manhattan, KS 66506. 



3t 



f^ ■ Department of Mechanical Engineering, Kansas State University, Manhattan, KS 66506. 

o 



Abstract 



o 

pi ^ A detailed analysis is presented to demonstrate the capabilities of the lattice Boltzmann 

S: 

O . method. Thorough comparisons with other numerical solutions for the two-dimensional, 

O 

driven cavity flow show that the lattice Boltzmann method gives accurate results over 
H , a wide range of Reynolds numbers. Studies of errors and convergence rates are carried 

out. Compressibility effects are quantified for different maximum velocities, and parameter 
ranges are found for stable simulations. The paper's objective is to stimulate further work 
using this relatively new approach for applied engineering problems in transport phenomena 
utilizing parallel computers. 



1 Introduction 

Lattice gas automata (LGA) and its later derivative, the lattice Boltzmann equation method 
(LBE), are relatively new approaches that utilize parallel computers to study transport phe- 
nomena. Since the first two-dimensional model representing incompressible Navier-Stokes 
equations was proposed by Frisch, Hasslacher, and Pomeau (FHP) in 1986 [1], LGA have 
attracted much attention as promising methods for solving a variety of partial differential 
equations and modeling physical phenomena [2, 3, 4, 5]. 

A lattice gas is constructed as a simplified, fictitious microworld in which space, time and 
velocities are all discrete. In general, a lattice gas consists of a regular lattice with particles 
residing on the nodes. A set of Boolean variables nj(x, t)(i = 1,- ■ ■ ,b) for describing the 
particle occupation is defined, where b is the number of directions of particle velocity at each 
node. Starting from an initial state, the configuration of particles at each time step evolves 
in two sequential steps: (a) streaming, where each particle moves to the nearest node in the 
direction of its velocity; and (b) colliding, which occurs when particles arriving at a node 
interact and possibly change their velocity directions according to scattering rules. For 
simplicity, the exclusion principle (no more than one particle is allowed at a given time and 
node moving in a given direction) is imposed for memory efficiency and leads to the Fermi- 
Dirac equilibrium distribution. The strategy of the lattice gas is two-fold: a) to construct a 
model as simple as possible of the microworld to permit simulations of a system composed 
of many particles and b) to capture the essential features of real collision processes between 
particles such that, for long times and large scales, macroscopic transport phenomena are 
captured. 

That the evolution of particles on an artificial lattice can simulate the macroscopic 



behavior of fluid flow is based on the foflowing facts: the macro-dynamics of a fluid is the 
result of the collective behavior of many particles in the system and details of the microscopic 
interactions are not essential. Changes in molecular interactions affect transport properties 
such as viscosity, but do not alter the basic form of the macroscopic equations as long as 
the basic conservation laws and necessary symmetries are satisfied [2, 3]. 

Due to the microscopic nature and local interaction between particles, the lattice gas 
approach possesses some unique advantages. The scheme is absolutely stable; boundary 
conditions are easy to implement; the model is ideal for massively parallel computing and 
the code is simple. The lattice gas method also contains some problems such as non- 
Galilean invariance, due to the existence of a density-dependent coefficient in the convection 
term of the Navier-Stokes equation, an unphysical velocity-dependent pressure and inherent 
statistical noise that requires a spatial (or time) averaging to obtain smooth macroscopic 
quantities. To avoid some of these inherent problems, several lattice Boltzmann (equation) 
models have been proposed [6, 7, 8, 9, 10, 11]. The main feature of the LBE is to replace 
the particle occupation variables, nj, (Boolean variables) by the single-particle distribution 
functions (real variables) /j = (n-j), where ( ) denotes a local ensemble average, in the 
evolution equation, i.e. the lattice Boltzmann equation. 

The lattice Boltzmann equation as a numerical scheme was first proposed by McNamara 
and Zanetti [6]. In their model, the collision operator is the same as in the LGA. Higuera, 
Jimenez and Succi [7, 8] introduce a linearized collision operator that is a matrix and has no 
correspondence to the detailed collision rules. Statistical noise is completely eliminated in 
both models; however, the other problems remain since the equilibrium distribution is still 
Fermi-Dirac. The lattice Boltzmann model proposed by Chen et al. [9, 11] and Qian et al. 



[10] abandons Fermi-Dirac statistics and provides the freedom required for the equilibrium 
distribution to satisfy isotropy, Gahlean invariance and to possess a velocity-independent 
pressure. Their models apply the single relaxation time approximation first introduced by 
Bhatnager, Gross and Krook in 1954 [12] to greatly simplify the collision operator. This 
model is called the lattice Boltzmann BGK model. 

Compared with the lattice gas approach, the lattice Boltzmann method is more compu- 
tationally efficient using current parallel computers. Applications have been done using both 
methods on hydrodynamics [13, 14, 15, 16], flow through porous media [17, 18], magneto- 
hydrodynamics [19, 20], multiphase flow [21, 22, 23, 24] and the reaction-diffusion equation 
[25, 26, 27]. Collected papers and applications of lattice gas and lattice Boltzmann methods 
can be found in [4, 5, 28, 29]. 

Despite these studies on various problems, thorough quantitative investigations of the 
method have not been published. In the present work, the lattice Boltzmann BGK model 
(LBGK) is used to solve for the viscous flow in a square, two-dimensional cavity driven 
by shear from a moving wall for Reynolds numbers up to 10,000. Detailed comparisons 
between the LBGK and traditional methods are presented. The compressibility error and 
the convergence rate of the method are discussed. The objective of this paper is to analyze 
the accuracy and physical fidelity of the lattice Boltzmann BGK method and to stimulate 
further studies using the lattice Boltzmann approach. 

Section 2 presents a technical synopsis of the lattice Boltzmann model used in this 
paper that will enhance the general reader's understanding of this simulation method. More 
technical details are given in the Appendix for those who want to use the lattice Boltzmann 
method. The lattice Boltzmann simulation of driven cavity flow is discussed in Section 3 



and thoroughly compared with results from other numerical methods. Section 4 studies the 
numerical errors in lattice Boltzmann simulations due to lattice size and compressibility. 
Section 5 is devoted to comparisons between the square lattice and the triangular (FHP) 
lattice. The limit of relaxation time for these two models is explored. The final section 
contains concluding remarks. 

2 Two-Dimensional Square Lattice Boltzmann Model 

In this section an outline is given of the procedures of the lattice Boltzmann simulation. 
A square lattice with unit spacing is used on which each node has 8 nearest neighbors 
connected by 8 links (see Figure lA in the Appendix). Particles can only reside on the 
nodes and move to their nearest neighbors along these links in the unit time step. Hence, 
there are two types of moving particles. Particles of type 1 move along the axes with speed 
|ei| = 1 and particles of type 2 move along the diagonal directions with speed |e2| = \/2. 
Rest particles with speed zero are also allowed at each node. The occupations of the three 
types of particles is represented by the single-particle distribution function, /o-i(x, t), where 
subscripts a and i indicate the type of particle and the velocity direction, respectively. The 
distribution function, f„i{x,t), is the probability of finding a particle at node x and time t 
with velocity eo-j. The particle distribution function satisfies the following lattice Boltzmann 
equation: 

fai{x + e„i,t + l) - f„i{x,t) =^ai, (1) 

where Q^i is the collision operator representing the rate of change of the particle distribution 
due to collisions. According to Bhatnagar, Gross and Krook (BGK) [12], the collision oper- 
ator is simplified by the single time relaxation approximation. Hence, the lattice Boltzmann 



BGK equation is 

^,(x + e,,,t + 1) - ^i(x,t) = -^[/..(x,t) - fif{^,t)], (2) 

where /^^ (x, f) is the equihbrium distribution at x, t and r is the single relaxation time 
which controls the rate of approach to equilibrium. The density per node, p, and the 
macroscopic velocity, u, are defined in terms of the particle distribution function by 

E E /- = p^ (3) 

17 i 

and 

E E J^'^i^'^i = P^- (^) 

o- i 

The equilibrium distribution can be chosen in the following form for particles of each type: 

#(0) 4 3 21 

/oi = 9^[^~ 2^ ^' 

/J? = ^P[l + 3(ei.-u) + ^(ei,.u)2-^n2], 
4? = ^^[1 + 3(e2. • u) + ^(e2. • u)2 - ^^^j, (5) 

The relaxation time is related to the viscosity by 

2r - 1 

" = ^r- «=' 

where u is the kinematic viscosity. The detailed derivation of the LBGK model is given in 
the Appendix. 

Having chosen the appropriate lattice size and the characteristic velocity for the LBE 
system, the viscosity, u, of the problem can be calculated for a given Re number and then 
the relaxation time is determined by the formula above. Starting from an initial state of 
/o-j(x, t), the density and velocity fields and hence the equilibrium distribution function 



can be obtained. In each time step, the updating of the particle distribution can be spht 
into two substeps: colUsion and streaming. It is irrelevant which one is the first for a long 
time run. The collision process at position x occurs according to the right hand side of 
the Boltzmann equation given as Eq. (0). The resulting particle distribution at x, which is 
the sum of the original distribution and the collision term, is then streamed to the nearest 
neighbor of x, x + Ba-i, according to the particle velocity eg-j. The updating procedure is 
terminated for steady state problems when certain criteria are reached. The method can 
be used for transient problems, but this will not be discussed in this paper. 

The boundary condition commonly used at the solid wall of a fluid simulation is the 
no-slip condition for which the velocities vanish at the wall. This is implemented in the 
lattice gas and lattice Boltzmann methods with the bounce-back rule in which all particles 
hitting the wall are reflected back in the direction from which they came. 

Another lattice model commonly used in two-dimensional lattice gas and lattice Boltz- 
mann simulations is the triangular lattice (FHP model) [1, 2]. This is a two-speed (0 and 1) 
model in which the lattice constant (link) is equal to one. Simulations of cavity flow are also 
performed in this paper using this model. Comparisons between FHP and square lattice 
are discussed in Section 5. Two commonly used models for three-dimensional simulations 
are the 24- velocity FCHC [2] and the 14- velocity cubic models [10, 14] . 

3 Cavity Simulation 

The problem considered is two-dimensional viscous flow in a cavity governed by the Navier- 
Stokes equations. An incompressible fluid is bounded in a square enclosure and the flow 
is driven by the uniform translation of the top boundary. The fluid motion generated in 



a cavity is an example of closed streamline problems that are of theoretical importance 
because they are part of broader field of steady, separated flows. The literature is abundant 
for this flow configuration that shows rich vortex phenomena at many scales depending on 
the Reynolds number, Re. Numerical methods for solving the Navier-Stokes equations are 
often tested and evaluated on cavity flows because of the complexity of the flows. 

Most numerical solutions of two-dimensional cavity flow [32-40] use a vorticity-stream 
function formulation and discretize the incompressible, steady linear or nonlinear Navier- 
Stokes equations by finite difference [32-35], multigrid [36, 38, 39] and finite element [40] 
methods and their variations [37]. Earlier work was reviewed by O. Burggraf [32] where his 
numerical solutions of the nonlinear Navier-Stokes equations for Reynolds number up to 400 
showed a large primary vortex and two secondary vortices in the lower corners. The later 
studies of A. S. Benjamin and V. E. Denny [35], U. Ghia, K. N. Ghia and C. T. Shin [36], 
R. Schreiber and H. B. Keller [37] show that tertiary vortices are formed near the bottom 
corners for higher Reynolds numbers. The present results using the lattice Boltzmann 
method are compared with those done by Vanka [38], Schreiber and Keller [37], Ghia et al. 
[36] and Zhou et al. [39]. Ghia et al. obtained numerical solutions up to Re=10,000 with 
a 257x257 grid using the coupled strongly implicit multigrid method and vorticity-stream 
function formulation. Their work is the most comprehensive study of cavity flow to date. 

The present simulation uses Cartesian coordinates with the origin located at lower left 
corner. The top wall is moving from left to right with velocity U. The cavity has 256 lattice 
units on each side. Initially the velocities at all nodes, except the top, are set to zero. 
The X- velocity of the top is U and the y-velocity is zero. Uniform initial particle density is 
imposed such that the moving particle 1 has a density fraction of d = § = 0.3 per direction. 



The moving particle 2 has a density fraction of j per direction, and the rest particle has a 
density of 4(i. Therefore, the total density per node is p = 2.7. Using the uniform density 
distribution and velocities given above, the equilibrium particle distribution function, /j, is 
calculated according to Eq. (^). The evolution of fi can then be found by a succession of 
streaming and collision-like processes. After streaming, the velocity of the top lid is reset 
to the uniform initial velocity. At the end of each streaming and collision process cycle, 
the particle distribution function, fi, at the top is set to the equilibrium state and the 
bounce-back boundary conditions are used on the three stationary walls. The two upper 
corners are singular points which are considered as part of the moving lid in the simulations, 
but tests shown there is little difference if these two points are treated as fixed wall points. 
The uniform velocity of the top wall used in the simulations is U=0.1. The compressibility 
effects are discussed in Section 4.4. The Reynolds number used in the cavity simulation 
is defined as Re =\]L/v, where U is the uniform velocity of the top plate, L is the edge 
length of the cavity and v is the kinematic viscosity that is related to the single relaxation 
time as given in Eq. (0). All the results are normalized to allow comparisons between the 
present work and other results based on a unit square cavity with unit velocity of the top 
boundary. 

Steady-state solutions for cavity flow are obtained using the lattice Boltzmann method 
for Re=10, 100, 200, 400, 1,000, 2,000, 5,000, and 7,500. The Re=10,000 case is also run 
on a 256x256 lattice, but steady state cannot be reached because bifurcation takes place 
somewhere between Re=7,500 and 10,000. The results for Re=10,000 oscillate between a 
series of different configurations. For this reason the results presented in this paper are 
those for Re up to 7,500. The dependent variables of stream function, velocity, pressure 



and vorticity are calculated using the particle distribution function, /j. The dependent 
parameter of the drag coefficient of the driving wall is discussed also. 

3.1 Stream function 

Figure 1 (a-g) shows plots of the stream function for the Reynolds numbers considered. It is 
apparent that the flow structure is in good agreement with the previous work of Benjamin 
and Denny [35], Schreiber and Keller [37] and Ghia et al. [36]. These plots give a clear 
picture of the overall flow pattern and the effect of Reynolds number on the structure of 
the steady recirculating eddies in the cavity. In addition to the primary, center vortex, a 
pair of counterrotating eddies of much smaller strength develop in the lower corners of the 
cavity at higher values of Re. At Re=2000, a third secondary vortex is seen in the upper 
left corner (it is generated at a critical Re of about 1,200 according to [35] and that agrees 
with the results of the present work). For Re >5,000, a tertiary vortex in the lower right 
hand corner appears. A series of eddies with exponentially decreasing strength in the lower 
corners has been predicted [36]. Due to the compressibility effect of the LBE (discussed in 
Section 4.4), the tertiary vortex shown in the Figure 1 (g) oscillates. 

For low Re (e.g. Re=10), the center of the primary vortex is located at the midwidth 
and at about one third of the cavity depth from the top (see Figure l.a). As Re increases 
(Re=100), the primary vortex center moves towards the right and becomes increasing cir- 
cular. Finally, this center moves down towards the geometric center of the cavity as the Re 
increases and becomes fixed in its x location for Re > 5,000. The movement of the vortex 
center location versus Re is shown in Figure 2 along with the results given by Ghia et al. 
[36]. 
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To quantify these results, the maximum stream function value for the primary vortex 
and the minimum values for the secondary vortices along with the x and y coordinates of 
the center of these vortices are listed in Table 1. All the results presented use a uniform 
top velocity U = 0.1 except for Re=100 where U=0.01 is used. The reason is discussed in 
Section 4.4. Also listed are results selected from previous work [36-39]. Previous results 
agree with each other for Re < 1000, but vary for higher values of Re. The results of the 
present work and that of Ghia et al. [36] for stream function values agree within 0.2 % for 
all values of Re (Re=2,000 data was not given by [36]). The locations of the vortex centers 
predicted by the lattice Boltzmann method also agree well with those given by Ghia et al. 
[36]. 

Unlike finite-difference or finite-element methods that start from the steady-state, partial 
differential equations, the present method is a direct modeling method that evolves into 
steady state. The time to reach steady-state depends on the lattice size, the values of Re 
and the driving velocity, U. For all the cases run in this paper, steady-state is reached 
when the difference between the maximum values of the stream function for successive 
10,000 steps (process cycles) is less than 10~^. Considering the kinetic, direct, compressible 
and unsteady nature of the lattice Boltzmann method, the excellent agreement with entirely 
different methods such as Ghia et al. [36] is quite encouraging. 

The minimum values of stream function and the center of the secondary vortex in the 
upper left corner for Re=5,000 and 7,500 are listed in Table 2. These results also show good 
agreement with Ghia et al. [36]. 
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3.2 Velocity profiles 

Velocity components along a vertical and horizontal center line for several values of Re are 
shown in Figure 3. The velocity profiles change from curved at lower values of Re to linear 
for higher Re values. The near linear profiles of the velocity in the central core of the cavity 
indicate the uniform vorticity region generated in the cavity at higher values of Re. These 
results agree with those from previous studies [34, 36, 37]. 

3.3 Vorticity 

The plots of vorticity in Figure 4 (a-g) show that the steady cavity flow within closed 
streamlines at high Re consists of a central, inviscid core of nearly constant vorticity with 
viscous effects confined to thin shear layers near the walls. Batchelor [30] predicted these 
results from his model for separated eddies in a steady flow. As the Re increases, several 
regions of high vorticity gradients (indicated by concentration and wiggle of the vorticity 
contours) appear within the cavity. The thinning of the wall boundary layers with increasing 
Re is evident from these plots, although the rate of this thinning is very slow for Re > 5,000. 
The values of vorticity at the center of the primary vortex for different Re are listed in Table 
3. These values closely agree with the results of Ghia et al. [36] and approach the analytical 
value of 1.886 for an infinite Reynolds number calculated by Burggraf [32] using Batchelor 's 
model. 

3.4 Pressure 

Figure 5 (a-g) displays the pressure deviation contours for the present simulations. Since 
only the pressure gradient appears in the Navier-Stokes equation, the values of pressure 
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can defFer linearly. These plots are in good agreement with the static pressure given by 
Burggraf [32] (note that the top wall in [32] moves from right to left which is opposite to 
that in the present simulation). The pressure in [32] is obtained by integrating the Navier- 
Stokes equation given the velocity field, while the pressure in the lattice Boltzmann method 
satisfies the equation of state of the isothermal gas where it is proportional to the density. 
The observed agreement between these very different approaches demonstrates that the 
lattice Boltzmann BGK model is valid for simulating incompressible flow. By examing the 
closed contours in the pressure plots, it is seen that the inviscid core grows with increasing 
values of Re. In the opposite limit of Re approaching zero, pressure becomes a harmonic 
function and the contours cannot be closed but must end on the boundaries [32]. The 
present results support this last requirement. 

3.5 Drag on the top 

The drag force and the drag coefficient of the moving wall is calculated here for Re values 
considered in the study. The stress on the moving wall is given by the Newton's formula as 

du 

where u is the x component of velocity and fi is the kinetic viscosity. The drag force on this 
surface, F^, is defined as 

Fd= I Tyxdx =, / a-— ax = a > Ax, 

7o "" Jo ^dy ^ ^ Ay 

where Ux is the grid number in the x direction, L is the length of the square cavity and 

Ay = Ax = , _-^^ are the spacing of the lattice. The drag coefficient is then written as 

Fd 



C, 
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where p is the average density and U is the velocity of the top. The drag coefficient 
decreases as Re increases, as found in other laminar flow configurations. This can be seen 
by introducing the dimensionless quantities 

, u , X , y 

The drag coefficient can then be expressed as 

Uj^fi^dx' _ fi ^ du' ^^, _ 1 [^du\^, 



'' ~ pU^L ~ Jo PUL dy' ~ Re Jo dy' ' 
The results of drag and drag coefficient for different values of Re are listed in Table 4 and 
the latter is plotted in Figure 6. The plot shows that the relation between drag coefficients 
and Re satisfies the above formula. There is no data available on drag and drag coefficient 
from other methods for comparison. 

4 Error Analysis 

4.1 Sources of errors 

There is no analytic solution for cavity flow. Results from the work described in this paper 
are compared with the numerical solutions obtained by several other methods. Differences 
are found between the results of previous work, especially for higher values of Re. Several 
of these authors state that the data for the secondary vortices are less reliable due to corner 
singularities and/or roundoff errors [37], mesh-size limitations [34] or because the values 
of the stream function in the corners are small and, in some cases, below the convergence 
accuracy of the calculations [38] . 

Before the final results listed in Table 1 were obtained, the results from lattice Boltzmann 
simulations were very close to the results given by Ghia et al. [36] for Re> 1,000 only. The 
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properties of the secondary vortices were less satisfactory for Re less than 1,000. The 
secondary vortex of the lower left corner for Re=100, whose stream function is a small 
quantity of the order of 10~^, was not detected by the lattice Boltzmann method used in 
the initial simulations. Also, the secondary vortex in the lower right corner for the same 
Re, whose stream function is of the order of 10~^, did not match corresponding results of 
other investigators. Although these are not major features, it was important to investigate 
the cause of these discrepancies. 

The theoretical assumptions of the present method are the Boltzmann transport equa- 
tion plus the single relaxation time approximation of the collision term. As long as the 
macroscopic properties of fluid vary slowly enough in space and time compared with mi- 
croscopic particle dynamics, collisions should maintain approximately the local equilibrium 
such that the assumptions of molecular chaos by Boltzmann and single relaxation time by 
BGK are valid for problems of fluid dynamics. The possible reasons for the "errors" in the 
present simulations may be categorized as follows: 

1. The simulations done for different values of Re on a 256 x 256 lattice used single 
precision arithmetic. It is possible that roundoff error could be accumulated. 

2. The small compressibility effect presented in the LBE simulations may cause differ- 
ences when compared with models where compressibility is zero. 

3. The lattice size used here may still be too coarse to resolve all the small scale 
phenomena. 

4. The time step at which the simulation is terminated may not be large enough to 
represent steady state. 

5. The integration methods used in the calculation of the stream function may introduce 
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errors. 

The errors caused by 1. and 4. can be avoided by using double precision floating-point 
arithmetic and running for longer times. These experiments did not change the result for 
weak vortices on the smaller scales. What follows are subsections investigating the rest of 
the error sources. 

4.2 Effect of lattice size 

To test the eff'ect of lattice size, simulations for Re= 1,000 are done on the following lattice 
configurations: 33 x 33, 65 x 65, 129 x 129, and 513 x 513. The driving velocity used 
is kept at f^ = 0.1. Two relative velocity errors were calculated according to the following 
formula: 



El 


Ex 


,y 1^1 - 


-uo\ 


+ 


\vi - 


■^^ol 




2—ix,y 


\uo\ 


+ 


\vo\ 


; 


F?.= 


\/S^ 


,y{ui - 


uo) 


2 + 


{vi 


-voY 



/Ex,,(^o)2 + {voy 

where u, v are the x and y components of the velocity, respectively. The subscript 0, 1 
indicate the 513 x 513 and the coarser grain lattice, respectively. Velocities on different 
grids are taken at corresponding positions, while the sums are taken over the entire lattice. 
El and E2 are two common relative velocity errors in the LI and L2 norm sense, respectively, 
based on the highest resolution lattice (513). 

The results for El and E2 are plotted logarithmically in Fig. 7 and listed in Table 5. 
The quantities El and E2 are calculated also according to the formula above where the 
and 1 indicate two successive lattice sizes. The results are similar to that shown in Fig. 7. 
It is clear from Fig. 7 that the convergence rate is approximately first order in space. This 
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result is different from other works [15, 41] where a second-order convergence rate is claimed. 
In [15] the problem is a decaying Taylor vortex flow with periodic boundary conditions (no 
solid wall) and in [41] exact boundary conditions for the particle distribution function are 
used instead of the bounce-back condition implemented at the wall in the present work. 
The first order convergence rate observed here may be due to the bounce-back condition 
used on the stationary walls. This observation is confirmed by a recent paper [42] in which 
the bounce-back boundary condition is shown to be a first-order approximation for a no-slip 
wall. It is shown also in [42] that higher order accuracy can be achieved by improving the 
implementation of boundary conditions. 

Better resolution is obtained as the number of lattice nodes increases. However, com- 
puter time grows with lattice number because more nodes are updated and the time to 
reach steady-state is much longer. The time steps required to reach steady state for differ- 
ent lattice sizes are listed in Table 6. Beyond these facts, a simulation run on a 513 x 513 
lattice for Re=100 did not improve the method's ability to predict the secondary vortex in 
the lower, left corner. 

4.3 Integration error 

The most important features of the cavity flow are related to the stream function. The 
stream function used by Ghia et al. [36] (and others) was the primary variable of the 
solution. In the lattice Boltzmann model, however, the primary variable is the particle 
distribution function, /j. The velocity at each site is calculated from fi and the stream 
function is obtained by integrating the velocities. 

To investigate the error caused by integration, three integration rules (rectangular. 
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trapezoidal and Simpson) are used for Re=100. The results from all three rules are of 
the same order of accuracy if the integrations are taken in the same direction (of course the 
trapezoidal and Simpson rules are higher-order integration schemes than the rectangular 
rule and hence produce less error). Significantly different results are obtained by integrating 
from the four different directions, namely integrate u along y from top to bottom, integrate 
u along y from bottom to top, integrate v along x from left to right and integrate v along x 
from right to left. Theoretically, they should all give the same value for the stream function. 
Prom a numerical point of view, the integration should be taken from the smaller scale, 
otherwise the smaller scale would be dawn into roundoff error. However, in this particular 
case, integration from bottom to top contains significant error. The reason can be seen from 
the error formula. If the trapezoidal rules is used: 

Hf) = [ fi^)dx « hllfia) + f{a + h) + --- + lf{b)] = /„(/), 



where h = — - and n is the number of subdivisions within [a,b], the integration error can 
be expressed as 

Enif) = I{f) - Inif) = -^^^^r(c„). 

Here c„ is a point between a and b. The above formula can be improved by the asymptotic 
error formula [43]: 

En{f)^-^[nb)-f\a)]. 

The factor, j^, in the present case is about 1.3 x 10^^. The error then depends on the 
derivatives at the end points. In the case of integration taken from top to bottom or from 
bottom to top, the two derivatives have opposite signs and the error is enhanced. In addi- 
tion, the value of the derivative on the top is large. Numerical tests show that integration 
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from top to bottom gives inaccurate values (1.96e~^ and 1.66e~^ for the minimum stream 
function on the lower left and right corners, respectively). The integration from the bottom 
to the top predicts another vortex of orderlO"^ on the upper left corner that does not exist 
for Re=100 and is caused by errors (due to both integration and compressibility). On the 
other hand, integration from left to right or from the right to left gives the same signs on the 
two end derivatives, therefore decreasing the integration error. Since the left corner vortex 
is smaller than that on the right, the integration taken from left to right gives better results 
(— 1.76e~^on left and — l.lOe"^ on right) than the integration taken from the right to the 
left (— 2.4e~^ on left and — 1.13e~^ on right). The conclusion is that using a trapezoidal 
rule and integrating the velocity component, v, from left to right gives the most accurate 
integral. The error is then of order of 10^^. The Simpson rule has about the same accuracy. 

4.4 Compressibility effect 

It has been shown that the present LBE model represents the Navier-Stokes equation in the 
incompressible limit (see Appendix). But in the LBE simulation, the density cannot be a 
constant (otherwise pressure change cannot be described). It is important to find the effect 
of compressibility on the present solution. 

One quantity that represents compressibility is the mean variation of density. The mean 
density is defined as 

where A^ is the total number of nodes. The mean variation of density is given by 



^^ijuiP-m 



N 
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For Re=100, this mean density fluctuation, A, is calculated for U = 0.1, U = 0.05, U = 0.01 
and listed in Table 7 along with the Mach number, M = ^, where c^ = 4= is the speed of 
sound for the present model. The table shows that 

A(t/ = 0.05) ^ J A {U = 0.1), 

and 

A(U = 0.01) « — A ([/ = 0.05). 
25 

These results agree with the known relationship [16] that A is proportional to M^. (This 

relation can be seen from the dimensionless incompressible Navier-Stokes equations.) 

The compressibility effect can also be examined for the cavity flow problem as follows. 

In the steady case, the continuity equation represented by LBE is 

V • (pu) = 0, 

due to a non-constant p. The velocity u does not satisfy the incompressible continuity 
condition given by 

V-u = 0. 

It is from this equation that the stream function can be defined using u = -^ and v = —-^, 
where il) is the stream function. There is actually no exact definition for the stream function 
in LBE. Given a discrete velocity field obtained from the LBE calculation, an approximation 
of the stream function for the incompressible flow with V • (u) = needs to be constructed. 
The stream function definition written as 

ip = —vdx + udy 
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is still used to calculate the stream function. When integrating in only the x- or y-direction, 

the integral becomes 

ry 
V' = / udy, 
Jo 

or 

rx 

ip = — vdx. 
Jo 

In the case of incompressible flow in a cavity, the boundaries coincide with the zero stream 
function. The integrals then take the form 

ip = udy = / vdx = 0, 



JO 

where L is the total length of the wall. For an incompressible model, integration of u (or 
v) along the y (or x) direction from one edge of the cavity gives a theoretical value of the 
stream function at another edge of zero. In the actual computations, the stream function 
at the wall will not exactly equal zero because of roundoff and integration error. Due to the 
additional effects of compressibility in the LBE method, if the stream function is calculated 
by integrating v from the left to the right edge of the cavity, the values of the stream function 
on the right wall would indicate the error caused by compressibility, roundoff and integration 
errors. Since the trapezoidal rule gives the same results for integrations taken from opposite 
directions if there is no roundoff error, the roundoff error is found by comparing the values 
of the stream function on the left and right wall taken from opposite directions (the other 
sources of error, compressibility and integration, are the same for these two integrals). This 
error is less than about 10^^. The integration error is of the order of 10^^ as discussed 
above. Therefore, the maximum and the mean value of '0 at the right wall can be computed 
as an indicator of error due to compressibility if this value is larger than 10^*^. The mean 
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and maximum stream function at the right edge of the cavity is defined, respectively, as 



g^^ n^Hn.,j)) 



and 

Sm = maxlV'ln^jj)!, 

3 

where Ux = Uy = 256 is the number of nodes in the x- and y-direction, respectively. These 
values are calculated for U = 0.1, U = 0.05 and U = 0.01 for Re=100 and listed in Table 8. 
Again, the results show that Sa and Sm are proportional to M^. The compressibility error 
calculated by this argument can be used as a quantitative measure of the compressibility of 
the LBE method. The change of compressibility error with Re is calculated for U = 0.1 and 
listed in Table 9. The error caused by compressibility does not vary much with Re number. 
Actually with increased Re the error is slightly decreased, but is still of the same order of 
magnitude. 

It is clear that the error caused by compressibility has about the same order effect as the 
small scale phenomena in the cavity flow for low values of Re. By choosing the direction of 
integration for the stream function carefully, predictions of the small vortices are obtained. 
Furthermore, the compressibility error can be reduced by using smaller velocities at the top 
wall. The results for Re=100 in Table 1 are calculated using U = 0.01, while other Reynolds 
number use U = 0.1. However, the time steps required to reach steady-state for smaller top 
velocity increases dramatically as seen in Table 10. To overcome the compressibility error 
in the present LBE model, a new incompressible LBE model for steady-state flow has been 
developed and will be published in another paper [44]. 
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5 Triangular Lattice (FHP) Versus Square Lattice 

Simulations for cavity flow are also carried out on a triangular lattice (FHP). There are two 
types of particles on each node of a FHP model: rest particles and moving particles with 
unit velocity ej along 6 directions. In analogy to the procedures used for the square lattice 
in the Appendix, the equilibrium distributions for the FHP model are given as 

-p(O) J 2 2 

Jo = do - pu = ap- pu , 
/f = d + \piie, ■ u) + 2(e, • u)^ - ^u'] = P^ + lp[(e, • u) + 2(e, • u)^ - ^u\ 
If the ratio of rest and moving particle is defined as 

the pressure is determined by the following isothermal equation of state: 

(l-a)/9 3 



and the speed of sound is 



9 1 — Q 



2 A + 6 

The viscosity is related to the relaxation time through an equation of the form 

2t- 1 

v = . 

8 

Theoretically, the relaxation time, r, cannot be lower than 0.5 for a positive viscosity. To 
reach higher Re, the relaxation time can be lowed. Tests on a 128 x 128 lattice with a 
maximum velocity of U = 0.1 show that a critical value for r exists. Above this value, 
the simulation is smooth and reasonable physical patterns for the cavity flow are seen in 
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the real-time plots. However, below this critical value, some nonphysical patterns appear. 
Further reduction in the value of r would cause the simulation to be terminated by numerical 
blow up. Define the critical value of r as the lowest limit for the relaxation time that gives 
physically correct results. This limit varies with the ratio, A, maximum velocity, U, and the 
problem studied. If A is increased, the speed of the sound will be decreased and the Mach 
number is then increased if the velocity is unchanged. Table 11 lists the lowest relaxation 
time, Tmim and hence the highest Re number, Rcmax, obtained for different A, along with 
their speed of sound, Cg, and Mach number, M. Table 11 shows that the highest Re can be 
increased by increasing A. However, the compressibility error is also increased. 

Table 12 lists results for the same conditions as in Table 11, but for slightly different 
initial boundary conditions on the density. In case 1 of Table 11, the initial density on each 
node of the wall is the same uniform distribution as that for interior nodes. In case 2 of 
Table 12, the initial density on the wall is set equal to zero. It is clear that when A = 1, 
the lowest r is much higher than that in Table 11, but for large A the differences of Tmin 
between these two cases are diminished. The rest particles in the LBE method play the 
role of a particle reservoir. When the macroscopic velocity is higher, rest particles can be 
turned into moving particle and vice versa. Higher values of A = do/d mean a larger fraction 
of rest particles in the density behave like a fluid that is less rigid and more flexible (the 
compressibility is higher). When A = 1 in case 2, the lowest value of r for stable results 
is 0.5668. However, setting r = 0.5667 would make the computation blow up immediately 
due to the large initial density gradient on the wall and the relatively small fraction of rest 
particles. The lowest limit of r does not depend on the lattice size. Changing the maximum 
velocity, U, does change the lowest r slightly. However, the highest Re numbers obtainable 
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by this approach are much lower than that for U = 0.1 (see Table 13 for case 2). 

The square lattice corresponds to A = ^ = 4. Tests on a 128 x 128 square lattice with 
the maximum velocity of f/ = 0.1 show that the value of r cannot be smaller than 0.507 
(Re=5485) for cavity flow. Hence a simulation run on a 256 x 256 lattice with U = 0.1 can 
reach Re=10,000 (r = 0.50768) which is about the highest limit of Re on this size of lattice 
for cavity flow. Using small U did not produce a further reduction of r. The square lattice 
is better than the FHP lattice in the cavity flow simulations because the former can reach 
higher values of Re than the latter for the same maximum velocity and lattice size. Since 
the boundaries of the cavity are fitted better using the square lattice than FHP lattice, the 
formation of the vortices is more gentle in the simulation process using a square lattice than 
a FHP lattice. The ranges of parameters presented in this section are consistent with the 
results of linear stability analysis of LBE method without boundaries [45] . 

6 Conclusions 

The lattice Boltzmann method is a derivative of the lattice gas automata method and 
therefore inherits from the LGA some of its advantages over traditional computational 
methods. It is parallel in nature due to the locality of the transport of particle information, 
so it is well suited to massively parallel computing. Due to the kinetic description of 
the lattice Boltzmann method, it is easy to handle the complex boundary conditions and 
properties of a fluid system, such as flow through porous media and multi-phase flow. One 
important improvement due to the LBE method is that it can fully recover the Navier-Stokes 
equations at the macroscopic level including Galilean invariance and a velocity-independent 
pressure. However, there is a trade-off. The lattice Boltzmann method no longer has the 
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pure Boolean operation and numerical stability guaranteed by LGA. 

Detailed study of the cavity flow problem using the lattice Boltzmann method has shown 
that the method is accurate compared with conventional methods using the same mesh size. 
This verification produces confidence to apply the method to other complex systems. All 
aspects of the present work such as boundary conditions, parameter ranges, lattice size and 
compressibility effects are important when the method is applied to other problems. The 
following remarks are in order: 

1. The proper implementation of the boundary conditions is crucial for the lattice 
Boltzmann simulation. Various boundary conditions such as periodic, particle bounce- 
back, wind tunnel and constant flux conditions are commonly used for different situations 
in LBE. It is important that the boundary conditions applied for the simulation represent 
the correct physical problem. In the cavity simulation, for example, besides the uniform 
top velocity and no-slip conditions on the wall, the mass must be conserved globally. Any 
violation of this restriction will produce nonphysical results. Be aware that some improper 
boundary conditions can give a qualitatively reasonable flow but lead to quantitatively 
incorrect results. 

2. The range of parameters for the model is explored for the cavity simulations. Param- 
eters such as the lattice size, maximum velocity, the ratio of rest and moving particle and 
the single relaxation time are adjustable in LBE. The lattice size should be chosen so that 
a good resolution for all scales in the problem can be obtained at an affordable cost. The 
maximum velocity used in a simulation should be properly small for a low Mach number 
and hence low compressibility requirement and for the validation of the equilibrium distri- 
bution which is an expansion of small velocity. For the Chapman-Enskog expansion to be 
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valid, the spatial gradients of density and velocity should be small also. Since the maximum 
velocity and lattice size are limited, the single relaxation time needs to be small to achieve 
the higher Reynolds numbers. It is found that the lowest relaxation time leading to stable 
simulations depends on the ratio of rest and moving particles, the maximum velocity and 
the problem. To obtain a reliable simulation, the relaxation time should be chosen not too 
close to the lowest limit for the problem under investigation. On the other hand, since the 
lattice spacing in the LBE is unit, the parameter, r, can be understood as a mean free path. 
Therefore, r should be small enough compared with macroscopic characteristic length scale. 
This is a necessary condition that the microscopic statistics of the LBE will approach the 
Navier-Stokes equations as shown in the multiscale expansion (see the Appendix). 

3. The compressibility effect may become important when physical quantities of the 
smallest scale in an incompressible flow is comparable to the compressibility error. Using 
a smaller maximum velocity can reduce this error. However, it is probably not practical to 
predict small scales on the order of 10~^ or smaller by the present LBE method. 

4. The square lattice is better than the triangular lattice (FHP) in two-dimensional 
simulations because the former can reach higher values of Re number for the same lattice 
size and maximum velocity. 

5. The computer time used in the simulations is not compared carefully with other 
methods since the LBE includes transient effects in this problem and hence is not economi- 
cal compared with the multigrid method. There is no doubt, however, that the method can 
simulate unsteady and other complex problems on a parallel computer with time compara- 
ble, if not superior, to other methods. 

Lattice gas and lattice Boltzmann methods are relatively new approaches for transport 
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phenomena. It is apparent that further research on both theoretical and practical aspects 
is needed. Implementation of higher order boundary conditions, models for better resolving 
the small scale phenomena, applications in new fields and to discover new physics, improve- 
ment of thermodynamical models and careful studies for three-dimensional geometries are 
challenges for future research. 
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Figure captions 

Figure 1. (a) Stream function for Re=10. Top velocity U=0.1 The center of the primary 
vortex is at (0.5216, 0.7686). The center of lower left vortex is at (0.0392, 0.0431). The 
center of lower right vortex is at (0.9647, 0.0392). 

Figure 1. (b) Stream function for Re=100. Top velocity U=0.01. The center of the 
primary vortex is at (0.6196, 0.7373). The center of lower left vortex is at (0.0392, 0.0353). 
The center of lower right vortex is at (0.9495, 0.0353). 

Figure 1. (c) Stream function for Re=400. Top velocity U=0.1. The center of the 
primary vortex is at (0.5608, 0.6078). The center of lower left vortex is at (0.0549, 0.0510). 
The center of lower right vortex is at (0.8902, 0.1255). 

Figure 1. (d) Stream function for Re=1000. Top velocity U=0.1. The center of the 
primary vortex is at (0.5333, 0.5647). The center of lower left vortex is at (0.0902, 0.0784). 
The center of lower right vortex is at (0.8667, 0.1137). 

Figure 1. (e) Stream function for Re=2000. Top velocity U=0.1. The center of the 
primary vortex is at (0.5255, 0.5490). The center of lower left vortex is at (0.0902, 0.1059). 
The center of lower right vortex is at (0.8471, 0.0980). 

Figure 1. (f) Stream function for Re=5000. Top velocity U=0.1. The center of the 
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primary vortex is at (0.5176, 0.5373). The center of lower left vortex is at (0.0784, 0.1373). 
The center of lower right vortex is at (0.8078, 0.0745). 

Figure 1. (g) Stream function for Re=7500. Top velocity U=0.1. The center of the 
primary vortex is at (0.5176, 0.5333). The center of lower left vortex is at (0.0706, 0.1529). 
The center of lower right vortex is at (0.7922, 0.0667). 

Figure 2. The locations of the center of the primary vortex for different values of Re 
numbers. The origin is the geometric center of the cavity. 

Figure 3. (a) Velocity profiles for u through the geometric center of the cavity. 

Figure 3. (b) Velocity profiles for v through the geometric center of the cavity. 

Figure 4. Vorticity contours of the cavity flow, (a) Re=10. (b) Re=100. (c) Re=400. 
(d) Re=1000. (e) Re=2000. (f) Re=5000. (g) Re=7500. 

Figure 5. Pressure deviation contours of the cavity flow, (a) Re=10. (b) Re=100. (c) 
Re=400. (d) Re=1000. (e) Re=2000. (f) Re=5000. (g) Re=7500. 

Figure 6. Drag coefficient of top wall versus. Re. 

Figure 7. Convergence rate of the LBE method for cavity flow at Re=1000 with top 
velocity U=0.1. The errors are calculated relative to results obtained on a 513 x 513 lattice. 
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Appendix 

This appendix details the derivation that shows how the Navier-Stokes equations are 
recovered from a lattice Boltzmann equation on a square lattice by using the Chapman- 
Enskog expansion procedure of kinetic theory. In addition, the equilibrium distribution 
functions are obtained to guarantee that the requirements of isotropy, Galilean-invariance 
and velocity-independent pressure are satisfied. 

On each node of a square lattice there are three types of particle, namely, a rest parti- 
cle, a particle moving along perpendicular directions and a moving particle along diagonal 
directions, (see Figure lA). 



Figure lA Schematic of a square lattice 



The velocity vectors ei^i,e2,j are defined as 



i — 1 i — 1 

ei j = (cos— — 7r,sm— — tt) i = !,■■■, 4, 

^2,1 = V 2(C0S(— -— TT -I- — ), Sin{ IT + —)) i = 1, • • • , 4. 
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The symmetric properties of the tensor 

i 

are needed in the derivation and given as follows: 
The odd order of tensors are equal to zero, i.e. 

^e,, = 0, a = 1,2, (1) 

i 

X! eaiaeail3eai-i = 0, «, /3, 7 = 1, 2, (2) 

i 

and 

X! eaiae-aipe-ai-ie-aiee-aic. = 0, a, /5, 7, 0, C = 1, 2. (3) 

i 

The second order of tensor satisfies 



X! e<Tiae<^i/3 = 2e^(5a/3, a, /3 = 1, 2, (4) 

i 



where 



3a,/3 



and 




is the length of eo-j . 

Finally, the fourth order tensor has an expression as 

4Aq,^^6I - SSa/S-yd, cr = 2, 
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where 

{1, a = (3 = "f = 0, 
0, otherwise, 
and 

Here the subscripts, a, denote the types of particle; i indicates the directions of particle 
movement and a, f3, 7, 6, and (" are the components of the coordinates. 

The Chapman-Enskog procedure is an asymptotic expansion method for solving the 
Boltzmann equation in kinetic theory. It is necessary to introduce a small parameter in an 
asymptotic expansion to compare orders. The lattice constants are required to be small 
compared with macroscopic characteristic length scale, for example, the edge length of the 
cavity, L, i.e. | eo-j | /L <^ 1. Using e as a measure of this small scale, instead of using a unit 
lattice constant and time step as in Section 2, the lattice Boltzmann equation is written as 

fai (x + ee^i ,t + e) - f„i(x,t) = Q^i , (6) 

where Q^i is the collision operator and e is solely used for distinguishing different orders. 
The lattice Boltzmann BGK equation is 

U (x + ee,, , t + e) - /,, (x, t) = - i [/,, (x, t) - /^f (x, t)] , (7) 

r 

where faiic = i=l; a = 1,2 i = l,---,4)is the single-particle distribution function, 
/^j (x, t) is the equilibrium distribution at x, t and r is the single relaxation time. A general 
form of /^j (x, t) can be taken as 

/i°^(x, t) = A^ + B^ie^i • u) + C^{e^^ • u)^ + D^u\ (8) 
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Here Aa-,Ba-,Ca and D^ are cofficients that need to be found and depend on p, but not 
on u. Equation (P) can be thought as a special type of small velocity (up to the u^ term) 
expansion of fl^ . 

It is obvious that Bq = Cq = and Eq. (^) can be written separately in the form 

fS\^,t) = Ai + B^ieu ■ u) + Ci{eu ■ u)^ + Diu", 

ff^{^,t) = A2 + Bi{e2i ■ u) + C2(e2. • u)^ + D^u", (9) 

where /^j and /o-i satisfy the following constraints: 



and 

.J 

<y i cr i 

These constraints mean that the non-equilibrium distributions do not contribute to the local 



EE^^ = EE/i? = P' (10) 



E E f<^^^^' = HY1 foi^^i = p^- (11) 



values of density and momentum. Using Eqs. (p^), ([Tl|), (|T|) and (^), some constraints for 
coefficients A^j, B^j, C^ and D^ are found to be 

Ao + 4^1 + 4^2 = P, (12) 

2Ci + 4C2 + Dq + ADi + 4D2 = 0, (13) 

and 

2Bi + 4^2 = P- (14) 

Starting from the LBE with BGK collision operator, the Navier-Stokes equations can 
be recovered. Taking a Taylor expansion of Eq. (j^ gives 

Uiix + ee„,,t + e)-Uii^,t) = ^ i^[- + (e,, • V)]"/,i(x, t), (15) 

ra=0 
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where 



eo-i • V — Ca-iada + e^ipOp 



and the Einstein summation convention is used. If terms up to O(e^) are retained in Eq (15) 
the results is 



6[| + (e.. • V)]/.. + ^[| + (e.. • V)fU + 0{e^) 



i[/.,(x,t)-/(?(x,t)]. 



(16) 



Next, the Chapman-Enskog-like expansion is applied to Eq.(|16[). Expanding /o-i about f^^ 



.(0) 



gives 



with constraints 



and 



u^ = Y.^C = a + ^f^ + ^f)^ + 



n=0 



HHf. 





' 
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.(0) 


1 
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at 










Gai 




pu 



The collision operator, ilo-i, becomes 



(n) 



0, n > 1. 



(17) 



(18) 



(19) 



^l^f!;^' + ^1!;^' + ■■■]■ 



(20) 



To discuss changes in different time scales, three time scales, to,ti,t2 are introduced as 



to — t, ti — et, t2 — e t, 



where 



d _ d d 2 d 

dt dto dt\ dt2' 



(21) 
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Substituting Eqs. ([l7|) , (po|), and (21) into Eq ([iq) , the equation to order of e is: 



(9,o+e.,.V)/i? = -i/«. (22) 

The equation to order of e^ is: 

5.o/iP + 5../i? + (e.. • V)/(p + i[|- + (e.. . V)]VS^ = -^/i?. (23) 

Using Eq. (^), one obtains 

(5to + e.. • V)Vi? = -\dto + ^.^ ■ V)/il^ (24) 



Substituting Eq. (pi) into Eq. (p3) leads to 



5*1 /i? + (1 - ^)(5to + e.. • V)/(l) = -i/i?. (25) 



To derive the equations for p and pu to the first order in e, a summation of Eq. ( p2D with 
respect to o" and i is taken to give 

9.0 E E /S^ + E E(e.. • v)/(? = -i E E /iP = 0. 

which is the first order continuity equation 

5ioP + V-(pu)=0. (26) 

Similarly, multiplying eg-i in Eq. ( [2^ ) and taking summation as above gives 



Qt. E E /i?^- + E E(^- • v)/i?e.. = -^ E E /i^- = 0' 

which can be simplified to 

9,„(pu)+V-^^(e.,e.,)/(? =0. (27) 

o" i 
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Defining the momentum flux tensor as 

n = EE(<^-<^-)/-' (28) 

cr i 

Eq. (27) can be rewritten as 



aio(pu) + v-n(°) = o. (29) 

Similarly, the equations of order of e^ for p and u can be obtained from Eq. (^) as 

9../^+(i-^MoEE/il^ + (i-^)v-EE/il^e.^ = -7EE/i? = 0' 



o- ?, 



and 



dt, ipu) + (1 - -)a„ E E /il^e- + (1 - ^)v • E E e-e-/il^ = -r E E fif^^^ = o- 



O" J 



Here the constraint given by Eq. ( |18[) is used. By applying another constraint, Eq. (|19[), 
these two equations can be simplified to 

dt,p = 0, (30) 



and 



at,(/,u) + (i-^)v-nW = o. (31) 

Write n*-'') as follows: 

a i 

— / ^ / ^ ^aia^ai(i\-^a + -D^jeuiyU-y + Uue.(ji^€Q-igU.yUg -\- -Uq^U j 
cr i 

= E ^TSe^^a/J + E CaU^ug E e(Tiae^i/3eai7ec.ie» + E ^a'^e^a^apu 

a a i a 

= {2Ai + 4.A2)5ap + Ciu^ue{25ai3ye) 
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Since 

and 

one can write n|^i as 

nJJ] = [2Ai + 4^2 + (4C2 + 2Di + AD2y]6af3 + SCsu^u^ + (2Ci - SCs)^^^/?^^/?. (32) 

The first term is the pressure term and the other two are nonhnear terms. In order to 
obtain velocity-independent pressure, the coefficient of u^ is chosen to satisfy 

4C2 + 2Di + 4D2 = 0. (33) 

To have Galilean invariance, the non-isotropic term is eliminated by choosing 

2Ci - 8C2 = 0. (34) 



Eq. (|32D becomes 



Assuming that 



nJJ) = {2Ai + 4A2)6^p + 8C2UaUp. (35) 



8C72 = p, (36) 



and 

2^1 + 4^2 = cIp, (37) 

where Cs is speed of sound, gives the final expression for n''^^ as 

^a^ = <^lP^al3 + PUaUp- (38) 
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Substituting Eq. (|3q ) into Eq. (^) results in 

^ + V.(puu) = -V(cf/,). (39) 

Eqs. (^) and (^) are Euler equations that are derived from the e-order of the expansion 
of the Boltzniann equation. 

To derive the equations accurate to e^, the quantity V-n'-^) needs to be evaluated. From 
Eq. (22) the non-equilibrium distribution can be expressed as 



fl]^ = -rdtjl':^-r{e,,,d,)fl'lK (40) 

Substituting /^^ into Il^i gives 

= -^{(^to^al + "^i X! X! ^(Tiaeaipeai-iiAa + B„e„ieue + Cae^ieeaisUeUs + D^u^)}. 
Using Eq. (^) for n)^i leads to 

n„j = -T{dt^,[{clp)5ap + pUaU/3] + d^BiUg25af^^g + d^B2Ug{4:Aa(3'ye " ^dafB-yo)} 
= -T{-C^Ja0d-y{pU^) + dtoipUaUis) + 9q(2Bi - 8B2)ui35af3 

+4d^iB2U^)6ap + 4dUB2Up) + 4dpiB2Uo,)}. (41) 

To avoid non-isotropy, set 

2Bi - 8B2 = 0. (42) 

Recalling Eq. (jlj), Bi and B2 can be uniquely determined as 

B2 = -^, Bi = ^. (43) 

12' 3 ^ ' 

Therefore, Eq. (^) can be written as 

n^^j = -T{-d^{pU^)5af3 + i:da{pU/3) + -df^ipUa) - cld^{pu^)5ap + OtoipUaUfs)} . (44) 
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The last term can be simplified using Eq. ( p9D to take the form 

dtoipUaUp) = UpdtoipUa) + pUadt^UfS 

= Ui3[-d^,{pUaU-y) - daic^p)] + UadtoipUp) - UaUpdi^^p 

= -Upd^{pUaU^) - Updaiclp) - Ua[d-y{pU/3U^) + <9/3(c^/9)] + UaU/sd^ipU^) 
= -Uadp{clp) - Ufjdaiclp) - Ua,Ufsd^{pU^) - pUaU-yd^Up + UaUpd-y^pU^) - Upd^{pUaU^) 
= -Uadi3{clp) - Updaiclp) - pUaU^d^Up - Ui3dj{pUaU:y). 



Eq. (44) therefore becomes 



n„^ = -r{-d^ipu^)Sap + -daipup) + -dpipua) - c,d./{pu^)6ap 

-Uadp{clp) - Updaic^p) - pUaU^d^Up - UpdjipUaUj)} . (45) 



Note that Il^i does not only depend on the first spatial derivatives of p and u. Ignoring 



the last two terms which are order of 0{u ) in Eq. (|45|), leads to 

+^df3{pua) - Uadpiclp) - updaiclp)} + 0{u^). (46) 

Combine equations of 0(e) and O(e^) for p and u, and Eqs. (|26|), (|39|), (|30| ) and (0T|) with 



Eq. (46) as follows: 



Eq. (|26D added to Eq. (|30|) multipled by e gives 

dtoP + edt,p + V ■ (pu) = 0, 

which gives the correct form of the continuity equation as 

dp 



^^+V-(pu)=0. (47) 
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Eq. (H) added to Eq. (|3l|) multipled by e gives 



dtipn) +V ■ (pun) = -V ■ (cip) - eil 



2t' 



)v-nS. 



(48) 



(1) 



Substituting Eq. (|4|) for U\J^, Eq. (||) becomes 



dtipua) + dfsipUaUfi) = -daic^p) + e(T - -)dp{{- - cl)dy{pu^)6ai3 
+-da{pui3) + -dp{pua) - u^dpiclp) - updaiclp)] + 0{v?) + O(e^), 



(49) 
(50) 



which may be written in the form 

1 1 

dt{pua) + dp{puaup) = -da{clp) + e(r - 2){^"[(3 " c^)a^(/)u^)] 

+dp[^p{d^up + 5;3W,)] + dp[{^ - cl){u^dpp + upd^p)]} + 0{u^) + 0(e=^). (51) 



Consider the constraints on A^ given in Eqs. (12) and (p7|), and choosing 



^0 = T,P, M = -p, A2 = —p, 



Eq. (12) is satisfied and the sound speed is determined by 



Eq. (51) is simphfied as 



.2„^ , 1.,, ^^„r..^..„ , ^^r.. M , ^..3^ 



dtipUa) + d(3{puaup) = -do^icip) + ^e(T - l)dp[p(d^up + dp{u^)] + 0{u') + 0(e3). (52) 



Define the stain-rate tensor as 



Sal3 = ^{daUfS + dpUc 



(53) 



Then Eq. (p^) can be rewritten as follows: 



dtipuc,) + dnipuc^Ufi) = -daicip) + 2i^dp{pSap) + 0{u^) + Oie^), 



(54) 
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where 

2r - 1 
V = —^6, (55) 

with V being the the kinematic viscosity. Recall the Navier-Stokes equations in the two- 
dimensional space [18] 

dtipUa) + dfsipUaUfs) = ~daP + 5/3{2/i(5a^ - -Ujy6afs)}, (56) 

and 

dtp + V- (pu) = 0. (57) 

For an incompressible fluid with constant viscosity, /i, the Navier-Stokes equations become 

dtipUa) + dpipUaUp) = -daP + df3{2nSa/3} = -daP + pdfsdpUa- (58) 

For p = constant, the incompressible Navier-Stokes equations are 

V-u = 0, 

dtUa + dpiUaUfs) = -5q,(-) + i^djua- (59) 



It is seen that Eq.(p4D is exactly the same as the Narier-Stokes equation ( pq ) in the incom- 
pressible limit, V • u = 0. 

Collecting all coefficients so far, one obtains 

^0 = gP, A, = -p, A, = -p, 

Bi = Ip, B, = Ip, 

The remaining coefficients Dq,Di and D2 are related by Eq. ( [T^ and Eq. p^), so there is 
one free parameter. Since all coefficients of particle 2 are one-fourth of the corresponding 
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coefficients of particle 1, one can require Di = AD2- Hence, the remaining coefficients are 
determined as 

Finially, tlie equilibrium distribution functions are given as 

^(0) ^ n ^ 2i 
/oi = 9^[^" 2'' ^' 



fii = nP[l + 3(ei* • u) + -(ei, • uf - -u% 



f(0)_ 1 M , Q/. ^ , 9^. A2 3^21 

/i? = ^P[l + 3(e2. • u) + ^(e2. • u)2 - ^u% (60) 
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